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Abstract 

Monte Carlo simulations of neutron-rich matter of relevance to the inner neutron-star crust are 
,_^ performed for a system of A = 5, 000 nucleons. To determine the proton fraction in the inner 

^~~i crust, numerical simulations are carried out for a variety of densities and proton fractions. We 

P^ conclude — as others have before us using different techniques — that the proton fraction in the 

,__( inner stellar crust is very small. Given that the purported "nuclear pasta" phase in stellar crusts 

, ^ develops as a consequence of the long-range Coulomb interaction among protons, we question 

_^ whether pasta formation is possible in such proton-poor environments. To answer this question, 

^-H we search for physical observables sensitive to the transition between spherical nuclei and exotic 

^ pasta structures. Of particular relevance is the static structure factor S{k) — an observable sensitive 

1 -Q to density fluctuations. However, no dramatic behavior was observed in S{k). We regard the 

,_I^ identification of physical observables sensitive to the existence — or lack-thereof — of a pasta phase 

in proton-poor environments as an open problem of critical importance. 
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I. INTRODUCTION 

Neutron stars are unique laboratories for the study of matter under extreme conditions 
of density and isospin asymmetry [H [2] . Indeed, the conditions in the interior of neutron 
stars are so extreme that they are unattainable in terrestrial laboratories. Thus, neutron 
stars — and the exotic phases within — owe their existence to the presence of enormous grav- 
itational fields. To maintain hydrostatic equilibrium throughout the star, these enormous 
gravitational fields must be balanced by the pressure support of its underlying constituents. 
This results in an enormous dynamic range of pressures and densities that enables one to 
probe the equation of state (EOS) far away from the equilibrium density of normal nuclei [3]. 

Neutron stars contain a non-uniform crust above the uniform liquid mantle (or stellar 
core). The core is comprised of a uniform assembly of neutrons, protons, electrons, and 
muons packed to densities that may exceed that of normal nuclei by up to an order of mag- 
nitude. The highest density attained in the core depends critically on the equation of state 
of neutron-rich matter which at those high densities is presently poorly constrained. The 
core accounts for almost all of the mass and most of the size of the neutron star. However, at 
densities of about half of normal nuclear density, the uniform core becomes unstable against 
cluster formation. At these "low" densities the average inter-nucleon separation increases 
to such an extent that it becomes energetically favorable for the system to segregate into 
regions of normal density (nuclear clusters) and regions of low density (dilute neutron va- 
por). Such a clustering instability signals the transition from the uniform liquid core to 
the non-uniform crust. Note, however, that the precise value of the crust-to-core transition 
density is presently unknown, as it is sensitive to the poorly constrained density dependence 
of the symmetry energy [5] . 

The solid crust is divided into an outer and an inner region. In particular, the outer crust 
spans a region of about seven orders of magnitude in density (from about lO^g/cm^ to 4 x 
lO^^g/cm^ [SHIl)- Structurally, the outer crust is comprised of a Coulomb lattice of neutron- 
rich nuclei embedded in a uniform electron gas. As the density increases — and given that the 
electronic Fermi energy increases rapidly with density — it becomes energetically favorable 
for electrons to capture into protons. This results in the formation of Coulomb crystals of 
progressively more neutron-rich nuclei. Eventually, the neutron-proton asymmetry becomes 
too large for the nuclei to absorb any more neutrons and the excess neutrons go into the 
formation of a dilute neutron vapor; this signals the transition from the outer to the inner 
crust. 

In this contribution we are interested in modeling the structure of the inner crust. The 
outer-to-inner crust transition density is predicted to occur at about 4 x lO^^g/cm^ [5H7]. At 
this — neutron-drip — density the neutron-rich nucleus (^^^Kr) that comprises the crystalline 
lattice is unable to retain any more neutrons. Thus, the top layer of the inner crust con- 
sists of a Coulomb crystal of neutron-rich nuclei immersed in a uniform electron gas and 
a dilute — likely superfluid — neutron vapor. In contrast, at the bottom layer of the inner 
crust the density has become high enough (of the order of lO^^g/cm^) to "melt" the crys- 
tal into a uniform mixture of neutrons, protons, and electrons. Yet the transition from the 
highly-ordered crystal to the uniform liquid is both interesting and complex. This is because 
distance scales that were well separated in both the crystalline phase (where the long-range 
Coulomb interaction dominates) and in the uniform phase (where the short-range strong 
interaction dominates) become comparable. This unique situation gives rise to "Coulomb 
frustration". Frustration, a phenomenon characterized by the existence of a very large num- 



ber of low-energy configurations, emerges from the impossibility to simultaneously minimize 
all "elementary" interactions in the system. Whereas protons are correlated at short dis- 
tances by attractive strong interactions, they are anti-correlated at large distances because 
of the Coulomb repulsion. Whenever these short and large distance scales are well separated 
(as in the outer crust) protons bind into nuclei that are then segregated in a crystal lattice. 
However, as these length scales become comparable — at densities of about 10^^—10^^ g/cm^ 
as in the inner crust — competition among the elementary interactions results in the forma- 
tion of complex topological structures collectively referred to as "nuclear pasta" . Given that 
these complex structures are very close in energy, it has been speculated that the transi- 
tion from the highly ordered crystal of spherical nuclei to the uniform phase must proceed 
through a series of changes in the dimensionality and topology of these structures [H [9]. 
Moreover, due to the preponderance of low-energy states, frustrated systems display an 
interesting and unique low-energy dynamics. 

Interestingly enough, a seemingly unrelated condensed-matter problem — the strongly- 
correlated electron gas — appears to be connected to the nuclear pasta. In the case of the 
electron gas, one aims to characterize the transition from the low-density Wigner crystal 
(where the Coulomb potential dominates) to the uniform high-density Fermi liquid (where 
the kinetic energy dominates) [lO] . It has been argued that instead of the expected first order 
phase transition, the transition is mediated by the emergence of "niicroemulsions" , namely, 
exotic ("pasta-like") structures. Indeed, it has been shown that in two-spatial dimensions 
first-order phase transitions are forbidden in the presence of long-range (e.g., Coulomb) 
forces [TT]. One should note that no generalization of this theorem to three dimensions 
exists. Hence, although the existence of pasta phases in both core-collapse supernovae and 
neutron stars may be plausible P, |9], there is no guarantee that they exist. Actually, in a 
recent publication Oyamatsu and lida have shown that pasta formation may not be universal 
and that its existence (or lack-thereof) is intimately related to the density dependence of the 
symmetry energy [12]. In particular, the authors concluded that pasta formation requires 
models with a soft symmetry energy, namely, ones that increase slowly with density. 

Although the complex dynamics of the electron gas may shed light on the possible emer- 
gence of pasta phases in the stellar crust, an essential difference between these two systems 
remains. Whereas the constituents of the electron gas (i.e., electrons) are all electrically 
charged and thus experience long-range forces, the neutrons in the crust interact solely via 
short-range forces. Thus, a large neutron fraction in the inner crust could hinder the for- 
mation of the nuclear pasta. Indeed, pure neutron matter does not cluster. So given that 
the proton fraction in the crust is highly sensitive to the density dependence of the sym- 
metry energy, the emergence of pasta phases involves a delicate interplay between: (i) the 
symmetry energy (which controls the proton fraction), (ii) the surface tension (which favors 
spherical nuclei), and (iii) the Coulomb interaction (which favors deformation). It is this 
interesting and unique problem that we propose to study here via numerical simulations. 

Monte-Carlo simulations will be performed directly in terms of the nucleon constituents. 
Although we have already performed simulations of this kind [T314T5] . in the present 
manuscript we aim to improve them in two critical areas. First, in our previous work a 
screened Coulomb interaction between the protons was assumed with a screening length 
fixed at 10 fni. Although it was argued that no major qualitative changes are expected from 
such an approximation, in this contribution we evaluate the Coulomb interaction exactly 
via an Ewald summation. For completeness, an appendix has been included where a de- 
tailed derivation of the Ewald summation is presented for a system of protons embedded 



in a uniform electron gas. Second, our earlier simulations were all performed at the single 
proton fraction of Yp = 0.2. This value was chosen as it is intermediate between the large 
values attained in core-collapse supernovae and the low proton fractions characteristic of the 
inner crust of cold — fully catalyzed — neutron stars. Instead, in this work we will simulate 
neutron-rich matter for a variety of densities and proton fractions, albeit for a relatively 
small number of nucleons {A = 5,000). In this way, the optimum proton fraction will be 
determined by imposing /3-equilibrium. 

We should note that for the large proton fractions characteristic of core-collapse super- 
novae {Yp = 0.3-0.5) there appears to be agreement that the transition from the ordered 
Coulomb crystal to the uniform phase must proceed via intermediate pasta phases (such 
as rods, slabs, tubes, etc.). This agreement has been reached whether one uses numerical 
simulations [T3HIH] or mean-field approximations with a variety of effective interactions 1191 - 
23]. What is unclear, however, is whether such exotic pasta shapes will still develop in the 
proton-poor environment of the inner stellar crust. For example, mean-field models that 
impose /3-equilibrium predict proton fractions at densities of relevance to the inner crust 
of only a few percent p!9] - [2Tl 123] . To our knowledge, no numerical simulations have been 
carried out with proton fractions less than Yp = 0.1 ^16j. Thus, our manuscript revolves 
around the following three fundamental questions: (a) Do numerical simulations support 
the very low proton fractions predicted by mean-field calculations? (b) Does the transition 
from a highly-ordered crystal of spherical nuclei to a uniform Fermi liquid must proceed via 
intermediate pasta phases even when the proton fraction is very small? (c) Can one identify 
a set of robust physical observables that are sensitive to the formation of the nuclear pasta? 

We have organized the manuscript as follows. In Sec. [TT] we introduce the model and 
develop the formalism necessary to carry out the Monte-Carlo simulations. Results are 



presented in Sec. |III| for the potential energy as a function of density and proton fraction. The 
optimal proton fraction is then determined by minimizing the total energy per nucleon of the 
system. The pair-correlation function (for both protons and neutrons) and the corresponding 
static-structure factor are computed with the aim of identifying sensitive observables to 
the formation of the nuclear pasta. Our conclusions and suggestions for future work are 



presented in Sec. IV 



II. FORMALISM 

We start this section by reviewing the semi-classical model that while simple, contains 
the essential physics of "Coulomb frustration", namely, competing interactions consisting of 
a short-range nuclear attraction and a long-range Coulomb repulsion [13]. We stress that 
the long-range Coulomb interaction among the protons represents the critical ingredient 
for pasta formation. Hence the importance of establishing whether the proton fraction in 
the inner stellar crust is large enough to help (or hinder) the formation of the exotic pasta 
structures. A variety of static observables will be computed using Monte-Carlo simulations 
having neutrons, protons, and electrons as the main constituents. As the system is simulated 
directly in terms of its basic constituents, there is no need to assume a-priori any specific 
nuclear shape, such as rods, slabs, tubes, etc. 

The total potential energy of the system includes both short-range nuclear and long-range 
Coulomb contributions. That is, 

V^(ri, . . . , r^) = VNuci(ri, . . . , r^) + Vboui(ri, . . . , r^) . (1) 



The nuclear potential is assumed to consist of a sum of "elementary" two-body interactions 
of the following form: 
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where the inter-particle distance has been denoted by r^ = |rj — rj| and the isospin of the 
nucleon as r — with r = l for a proton and r = — 1 for a neutron. The model parameters (a, 
6, c, and A) were calibrated in Ref. y^ and are given by the following values: a = 110 MeV, 
h = —26 MeV, c = 24 MeV, and A = 1.25 fm^. Although not accurately calibrated, this 
set of parameters were fitted to several ground-state properties [I3] . In particular, Monte 
Carlo simulations using this simple interaction are able to account for the saturation of 
symmetric nuclear matter while precluding the binding of pure neutron matter at all densi- 
ties [13]. Note, however, that given that such classical simulations are unable to reproduce 
the momentum distribution of a quantum Fermi gas, the effective A^A^ interaction must 
be properly adjusted to compensate for this shortcoming. In particular — unlike a realistic 
neutron-neutron [nn) interaction that is attractive at intermediate distances — the effective 
nn interaction adopted here is (essentially) repulsive at all distances; see Fig. [I] 




FIG. 1: (Color online) Nucleon-nucleon potential employed in our semi-classical simulations |13] . 



For the Coulomb interaction we assume that the protons are immersed in a uniform 
neutralizing electron background. Note that at the densities of relevance to the stellar crust, 
the electrons behave to a good approximation as a non-interacting free Fermi gas. For such a 
system the Coulomb energy may be written in terms of the electrostatic potential as follows: 
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Here Z is the number of protons (and neutralizing electrons), V is the simulation volume, 
and the electrostatic potential $(r) satisfies Poisson's equation. Given that numerical simu- 
lations must involve a finite number of particles, we rely on periodic boundary conditions in 
an attempt to minimize finite-size effects. For the case of the short-range nuclear interaction, 
the implementation of the minimum image convention is both simple and accurate: a given 
particle in the system interacts only with the closest image of all remaining particles [231 - 126] . 
Indeed, the short-range character of the A^A^ force guarantees that the interactions with all 
images that are farther away will be exponentially suppressed. This prescription, however, 
is not valid for the long-range Coulomb potential. In this case a given particle in V interacts 
not only will all remaining particles in the simulation volume but, in addition, with all pe- 
riodic images. As such, the source of the electrostatic potential $(r) in Poisson's equation 
consists of all protons in the simulation volume V, all their periodic images, and the uniform 
electron background. That is. 



j=l n 



where n = {nx,ny,nz) is a triplet of integers and L^V^'^. 

A method that computes accurately and efficiently the Coulomb potential is based on the 
(90 year old!) Ewald summation [27j. The basic idea behind the Ewald sum is to add and 
subtract Z individual smeared charges at the exact location of each (point) proton. The role 
of each negative charge is to screen the corresponding (point) proton charge over a distance 
of the order of the smearing parameter a. As long as a is significantly smaller than the box 
length L, the resulting (screened) two-body potential will become short ranged and thus 
amenable to be treated using the minimum-image convention. What remains then is a peri- 
odic system of smeared positive charges together with the neutralizing electron background. 
In configuration space this remaining long-range contribution is slowly convergent. The 
great merit of the Ewald sum is that the long-range contribution can be made to converge 
rapidly if evaluated in momentum space, i.e., as a Fourier sum. Indeed, the Fourier sum 
is rapidly convergent because momentum (k) modes satisfying |k|a^ 1 make a negligible 
contribution to the sum. Hence, by suitably tuning the value of the smearing parameter a, 
the evaluation of the Coulomb potential may be written in terms of two rapidly convergent 
sums; one in configuration space and one in momentum space. The derivation and imple- 
mentation of the Ewald sum is now part of the standard literature [21] - [2Bt [25] . Yet, for both 
convenience and completeness, we have provided in the appendix a detailed derivation of the 
Ewald formula for a system of Z protons immersed in a neutralizing and uniform electron 
background. Here we only summarize the essential results. 

Assuming Z protons confined to a simulation volume V = L^ and immersed in a uniform 
neutralizing electron background, the Coulomb potential may be written in the following 
form: 

Vboui(ri, ...,rz)= i—j UcoAsi, ...,Sz). (5) 

This structure reveals that the Coulomb potential is an interaction with no intrinsic scale. 
That is, once a dimensionful parameter has been identified — vo = e'^/L in the present case — 
the dimensionless Coulomb potential U depends exclusively on the scaled coordinates Sj = 
Ti/L. As mentioned above, the Ewald construction relies on the introduction of an artificial 
smearing parameter a that naturally divides the Coulomb potential into two contributions: 



a short-range contribution that converges rapidly in configuration space and a long-range 
one that converges rapidly in momentum space. That is, 

t^Coul(Sl, . . . , S^) = - ^ iMsr(Sij) + Mlr(Sij) + Uq . (6) 

Here the — two-body — short- and long-range contributions are given by 

erfc(g/go) 
Msr(s) = , (7a) 
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where so = a/L is the dimensionless smearing parameter and Uq is an overall constant: 
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The short-range contribution [Eq. (7a)] represents a (gaussianly) screened two-body 
Coulomb potential. Indeed, the characteristic 1/s fall-off of the Coulomb potential is mod- 
ified by the complement of the error function which falls rapidly to zero for inter-particle 
separations significantly larger than the screening length sq. The long-range contribution to 
the Coulomb potential is customarily written in momentum space as a sum that involves 
the square of the charge form factor times a suitably modified Coulomb propagator in mo- 



mentum space [see Eq. (A. 22)]. Rapid convergence of the momentum sum is also achieved 



through the introduction of the smearing parameter [see Eq. (7b)]. However, the evaluation 
of the momentum sum often remains the most numerically demanding part of the compu- 
tation, so various efficient methods have been designed for the purpose of expediting the 
Fourier sum |25]. In our case we will settle for a relatively simple approach that ultimately 
leads to the form of uir(s) given above (see appendix). Not surprisingly, this form also in- 
volves a numerically demanding Fourier sum. Yet, the advantage of the present method is 
that Uir(s) may be precomputed in a fine mesh before the start of the simulation. Hence, 
during the actual simulation the evaluation of Mir(s) is implemented via a look-up table and 
a simple interpolation scheme [28j. This approach results in significant savings in processing 
time with little compromise in accuracy. In the present contribution we have selected a 
smearing parameter that is significantly smaller than the box length, namely, so = 0.12. 

III. RESULTS 

In this section we present results for the potential energy per nucleon as a function of 
density and proton fraction as obtained from our Monte Carlo simulations. In an effort to 
simulate the quantum-mechanical zero-point motion of the particles, the simulations were 
carried out at a fixed temperature of T = 1 MeV [T314T5] . Moreover, in all simulations the 
number of baryons was fixed at A = 5, 000 with the proton fraction varying from 1^ = to 
1^ = 0.25; this is contrast to our earlier simulations that assumed a fixed value of 1^ = 0.2. 
Having fixed the number of baryons, the proton fraction, the temperature, and the density 
of the system, a Metropolis algorithm was used to generate a total of 3,000 Monte Carlo 



sweeps — with each sweep consisting oi A = 5, 000 particle moves. The first 2,500 sweeps were 
used to therniahze the system and the last 500 to compute various physical observables. This 
was done "off-line" by first saving the coordinates of all 5, 000 particles to file (for each of 
the final 500 sweeps) and then using a post-processor to compute the observables. We focus 
initially on the energy per particle as a function of density and proton fraction as we aim to 
estimate whether the proton fraction in the bottom layers of the inner stellar crust is small 
enough to hinder the formation of the pasta. 
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FIG. 2: (Color online) Monte-Carlo results for the total (nucIear-plus-Coulomb) potential energy 
of a system of A = 5, 000 nucleons as a function of density and proton fraction. 

The potential energy per nucleon {V)/A as a function of density and proton fraction 
y is displayed in Fig. [2] (Note that we have used the symbol "?/" to denote an arbitrary 
value for the proton fraction; we reserve the symbol "1^" to denote the value obtained from 
applying the condition of /3-equilibrium) . The strong-nuclear interaction favors symmetric 
N = Z (i.e., y = l/2) clusters so (V^) is a monotonically decreasing function of y (at least until 
?/~ 1/2). The y = values are related to the equation of state of pure neutron matter and 
display the density dependence predicted by the present model. To enforce /3-equilibrium it 
is convenient — as well as accurate — to fit the y dependence of the potential to a quadratic 
form. That is, 

{V)ip,y) 



A 



Mp) + Mp)y + Mp)y' 



(9) 



The best-fit values obtained for the three coefficients are displayed in Table [Tj 



To determine the optimal proton fraction we compute the total energy per nucleon of 
the system {E/A) as a function of both density and proton fraction y. The optimal proton 
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TABLE I: Quadratic fit to the y-dependence of {V)/A as given in Eq. ^. 



fraction (1^) is then obtained by minimizing E/A with respect to y at fixed density. We 
write the energy per nucleon of the system in the following form: 
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Tp, and 



Here Am = mn — 'mp — me is the neutron-proton-electron mass difference and T„, 
Te represent the kinetic energy of the three constituents. In the above equation the ^'ap- 
proxiniate" sign has been used to indicate that the contribution from Am is negligible (see 
Fig. |4]) and that the electronic kinetic energy Tg — assumed to be that of a relativistic elec- 
tron gas — is much larger than the classical proton contribution {3yT/2). Note that the 
energy {kinetic plus rest mass) per particle of a relativistic Fermi gas of particles of mass 
m, and number density p = kp/Sir"^ may be expressed in terms of the dimensionless Fermi 
momentum XF = kp/m, (and yp = ^yl + x^) in the following compact form: 
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(11) 



Unfortunately, the neutron contribution T„ to the energy of the system is difficult to esti- 
mate. This is because at the densities of relevance to the inner crust, the neutrons are either 
part of the heavy clusters or part of the dilute neutron vapor. (See Fig. [3] for a snapshot of 
a Monte Carlo simulation that clearly illustrates this scenario.) Presumably, the neutrons 
bound to the heavy clusters are "sluggish" and thus may be treated classically. The "free" 
neutrons on the other hand, should probably be modeled as a quantum Fermi gas. However, 
deciding what fraction of the neutrons are bound to clusters and what fraction remains in 
the vapor is clearly a model-dependent question. Still, we can make some progress in setting 
reasonable limits for the proton fraction because the electronic Fermi energy dominates the 
kinetic energy of the system. To set a lower limit we assume no contribution to the kinetic 
energy from the neutron vapor. Conversely, a suitable upper limit for the proton fraction is 
obtained by assuming that all neutrons contribute to the kinetic energy. 

We illustrate how this process is implemented in Fig. |4] where the various contributions 
to the energy per nucleon have been plotted as a function of the proton fraction for a fixed 
density of p = 0.05 fm~^. As alluded earlier, the contribution from the mass difference Am 




FIG. 3: (Color online) Two snapshots from a Monte Carlo simulation for a system of yl = 4, 000 
nucleons at a baryon density of /? = 0.01 fm~ (upper panel) and /) = 0.025 fm~ (lower panel). In 
both cases a proton fraction of Z/A = 0.2 and a temperature of r = l MeV were used [3]. 



is negligible. Hence, in the limit of no contribution from the neutron vapor, the proton 
fraction emerges from a competition between the potential energy — which favors symmetric 
matter — and the electronic Fermi energy — which favors pure neutron matter. Given that the 
electronic Fermi energy changes more rapidly with proton fraction than (V), /3-equilibrium 
is reached for the very small value of yp~ 0.006. As expected, adding the vapor contribution 
(as a free Fermi gas of A^ neutrons) increases the proton fraction by almost a factor of 4; to 
y^~ 0.025. Still, in both cases the proton fraction is very small and significantly smaller — by 
almost a factor of 10 — than the 0.2 value employed in our earlier simulations [T3H15]. 

Having implemented the above procedure at all densities, we display in Fig. [5] the pre- 
dicted proton fraction as a function of baryon density. Our results suggest very low proton 
fractions over the entire range of densities explored in this work. This finding appears con- 
sistent with predictions made using vastly different approaches p^]42T| 123] . For example, 
in Fig. Is] results are also shown using the relativistic mean-field (RMF) model of Ref . |23] . 
The precipitous decline in the proton fraction displayed in the figure is driven by the rapid 
increase with density of the electronic contribution. That is, at the densities of relevance 
to the inner crust, electron capture is favored because the dilute neutron vapor contributes 
little to the energy whereas the electrons contribute a lot. This represents robust physics 
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FIG. 4: (Color online) Enforcing beta equilibrium at a density of p 
contributions to the energy per nucleon are explained in the text. 
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that should be fairly model independent. Thus, we must conclude that the proton fraction 
in /3-stable, neutron-star matter must be very small at the densities of relevance to pasta 
formation. Motivated by this finding, we find imperative to determine whether pasta for- 
mation can occur in the proton-poor environment of the inner stellar crust. To do so, one 
must be able to identify observables that are particularly sensitive to the formation of these 
exotic structures. For reasons that will soon become clear, we focus on the static structure 
factor. 

A fundamental physical observable that contains all information about the excitations 
of the many-body system is the dynamic response function S'(k, w). The dynamic response 
function — which is proportional to the scattering cross section — represents the probability 
that the system be excited by a probe (e.g., electron, neutrino, etc^j that transfers momen- 
tum k and energy u into the system. To obtain dynamical information of this kind one 
most resort to Molecular Dynamics simulations — which we hope to carry out in the near 
future. Unfortunately, Monte Carlo simulations such as the ones implemented here can only 
reveal static properties. Yet the static structure factor S{k.) — obtained by integrating the 
dynamic response over uj — provides a particularly useful (static) observable that is associ- 
ated with the mean-square density fluctuations in the ground state ^0\. Moreover, S{k) is 
intimately related to a quantity that can be readily determined in computer simulations: 
the pair-correlation function g{r). Indeed, ^(k) and g{r) are simply Fourier transforms of 
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FIG. 5: (Color online) Proton fraction computed under the assumption that either all the neutrons 
are bound into heavy clusters (iVyapor = 0) or all the neutrons contribute to the Fermi gas energy 
(A^vapor = -^) (see text for details). Also shown is a relativistic mean-field prediction for the proton 
fraction as described in Ref. 1231. 



each other. That is 



S{k) = l + yld\(^g{r)-l 



^— ikr 



(12) 



The pair-correlation function g{r) represents the probability of finding a pair of particles 
separated by a distance r. For a uniform fluid containing A^ particles and conflned to a 
simulation volume V, g{r) (which is only a function of the magnitude of r) may be computed 
exclusively in terms of the instantaneous positions of the particles. That is, 
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where the ^^brackets" represent an ensemble average. Whereas for a uniform fluid the one- 
body density is constant, interesting two-body correlations emerge as a consequence of the 
inter-particle dynamics and/or quantum statistics. For example, the characteristic short- 
range repulsion of the A^A^ interaction precludes particles from approaching each other. This 
results in a pair-correlation function that vanishes at short separations (see Fig. [6]). 

Given that the static structure factor accounts for the mean-square density fluctuations in 
the ground state, it becomes a particularly useful indicator of the critical behavior associated 
with phase transitions — which themselves are characterized by the development of large (i.e., 
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macroscopic) fluctuations. Indeed, the spectacular phenomenon of "critical opalescence" in 
fluids is the macroscopic manifestation of abnormally large density fluctuations — and thus 
abnormally large light scattering — near a phase transition [29j. In this regard, the static 
structure factor at zero-monientuni transfer provides a unique connection to the thermody- 
namics of the system [22]. That is, 

where ks is Boltzmann's constant, fi is the chemical potential, and k^ is the isothermal 
compressibility of the system. The isothermal compressibility kt is reminiscent of the spe- 
cific heat which accounts for energy, rather than density, fluctuations. They both play a 
fundamental role in identifying the onset of critical phenomena. 

We conclude this section by displaying in Fig. [6] the pair correlation function g{r) (upper 
panels) and static structure factor S{k) (lower panels) for both protons (left panels) and 
neutrons (right panels) at the fixed density of p = 0.03 fm~'^. At this value of the density 
pasta structures are known to be present — at least for moderately large proton fractions 
(see Fig. [3]). As expected, both (7„(r) and gp{r) vanish at small separations due to the 
characteristic short-range repulsion of the A^A^ interaction. Also in both cases, the first 
peak is particularly large given that nearest neighbors are strongly correlated within each 
cluster (see Fig. [s]). However, we also see qualitative differences between them. In the 
case of the neutrons, gn{r) displays the characteristic oscillatory structure of a fluid that 
is associated with the presence of nearest neighbors, next-to-nearest neighbors and so on. 
Unlike gn{f)-, however, the proton correlation function gp{r) displays only two prominent 
peaks — one large and one small. This is due to the absence of a proton vapor, as all protons 
are contained within neutron-rich clusters. Whereas the first peak in gp{r) develops due 
to the strong correlation between protons in the same cluster, the second (significantly 
smaller) peak represents the presence of protons in the neighboring clusters. This second 
peak appears to grow appreciably with proton fraction but then saturates. Yet it is unclear 
at this time whether the initial increase and eventual saturation may be significant. We will 
continue to explore this issue in the near future. 

The two lower panels in Fig. |6] display the static structure factors obtained as the Fourier 
transform of the corresponding pair-correlation functions. The static structure factor for the 
case of pure neutron matter {Yp = 0) is shown for reference and remains essentially struc- 
tureless for all values of k. In contrast, S{k) for neutron-rich matter displays a prominent 
peak that becomes progressively higher with increasing proton fraction. The maximum in 
S{k) reflects the most prominent oscillatory structure in g{r) which, in turn, is associated 
with the spatial (two-body) correlations in the system. Alternatively, the maximum in S{k) 
occurs at that momentum transfer for which the probe (e.g., electrons) can most efficiently 
scatter from the density fluctuations in the system. As in the case of the second maximum 
in gp{r), here too the maximum in Sp{k) grows with proton fraction and then appears to 
saturate. Whether this is significant in characterizing the emergence of the nuclear pasta 
remains to be investigated. What is clear, however, is that the static structure factor at 
zero momentum transfer S{k = 0) does not show the required enhancement associated with 
the onset of a phase transition. 
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FIG. 6: (Color online) Pair correlation function g(r) and static structure factor S{k) for protons 
(left panels) and neutrons (right panels) for a variety of proton fractions at a fixed baryon density 
of p = 0.03fm~3. 

IV. CONCLUSIONS 

Uniform neutron-rich matter at sub-saturation density is unstable against cluster for- 
mation. At densities below the neutron-drip line the ground-state of neutron-rich matter 
consists of a Coulomb crystal of spherical (neutron-rich) nuclei embedded in a uniform sea 
of electrons. It has been speculated that the transition from the highly-ordered crystal to 
the uniform Fermi liquid is mediated by an exotic phase known as "nuclear pasta". The 
nuclear-pasta phase — believed to exist in the inner crust of neutron stars — is characterized 
by the emergence of exotic nuclear shapes of various topologies that coexist with a dilute 
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vapor of neutrons and electrons. Critical to the formation of the pasta is the long-range 
Coulomb interaction among protons, as it promotes the deformation of the clusters and the 
development of exotic topological structures (see Fig. |3]). Given this essential fact, it is 
important to establish whether the proton fraction in the inner stellar crust is large enough 
to favor pasta formation. 

To compute the proton fraction in the inner stellar crust we relied on Monte Carlo simu- 
lations assuming that nucleons interact via a two-body, short-range nuclear interaction and 
a long-range Coulomb repulsion. We employed a simple A^A^ interaction identical to the one 
introduced in Ref. [13] that was (approximately) fitted to a few ground-state properties of 
finite nuclei, the saturation properties of symmetric nuclear matter, and that precludes the 
binding of pure neutron matter. Yet we improved on the work of Ref. [13] by treating the 
long-range Coulomb interaction exactly via an Ewald summation. Monte Carlo simulations 
were performed for a system oi A = 5, 000 nucleons at a variety of densities and proton 
fractions (note that the simulations in Ref. [I3] were limited to the single proton fraction 
of 1^ = 0.2). By doing so, the optimum proton fraction (at each density) was determined 
by imposing /3-equilibrium. The optimum proton fraction emerges from a delicate competi- 
tion between the nuclear symmetry energy — which favors proton fraction oi Yp< 1/2 — and 
the Fermi energy of the electrons — which favors Yp = 0. For the model employed here, the 
rapid increase of the electronic contribution with proton fraction dominates over the nuclear 
contribution. This leads to very small proton fractions at all densities of relevance to the 
inner crust. For example, at a density of p = 0.03 fm~^ the proton fraction is yp<0.02 — at 
least a factor of 10 smaller than assumed in Ref. [13]. Such small values for the proton frac- 
tion are consistent with those obtained using vastly different approaches (see, for example, 
Refs. [T9Vl2ni23] and references therein). Ultimately, the small proton fractions obtained in 
most approaches appears to be a direct consequence of the rapid increase with density of 
the electronic contribution. 

Given these results and the predominant role played by the Coulomb interaction in the 
formation of the nuclear pasta, it is only natural to ask whether such exotic structures 
can develop in the proton-poor environment of the inner stellar crust. To answer this 
question it becomes critical to identify observables that may be sensitive to pasta formation. 
Particularly useful in this case are observables sensitive to fluctuations, as these tend to 
increase rapidly near phase transitions. In this contribution we have focused on the static 
structure factor S{k) which is sensitive to density fluctuations. Moreover, the static structure 
factor is interesting as it provides a natural meeting place for theory, experiment, and 
computer simulations. Theoretically, S{k) is the Fourier transform of the pair-correlation 
function g{r) — a quantity that is amenable to computer simulations. Experimentally, the 
proton (neutron) static structure factor can actually be measured in electron (neutrino) 
scattering experiments. In this work we have computed the pair-correlation function — for 
both protons and neutrons — as a function of proton fraction for a fixed density of p = 
0.03 fm~ . In turn, the static structure factors were obtained by performing the required 
Fourier transforms. Although both g{r) and S{k) display interesting behavior that may be 
indicative of significant structural changes in the system, we find no clear evidence either in 
favor or against the formation of the nuclear pasta. In particular, the static structure factor 
at zero momentum transfer S{k = 0) does not display the large enhancement characteristic 
of a phase transition. Perhaps what is required to clearly identify the onset and evolution of 
the pasta phase are dynamical observables, such as various transport coefficients. To do so, 
one must rely on Molecular Dynamics (rather than Monte Carlo) simulations. A calculation 
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along these lines was reported in Ref. [30] where both the shear viscosity and thermal 
conductivity were computed. Perhaps surprisingly, no dramatic increase in the viscosity 
was observed as the system evolves from spherical to exotic nuclear shapes. Clearly, more 
effort should — and will — be devoted along these lines. 

In conclusion, our work revolved around three fundamental questions. First, what are 
typical values for the proton fraction in the inner stellar crust? Second, can exotic pasta 
structures develop in a proton-deficient environment? Finally, which physical observables 
are particularly sensitive to pasta formation? We have found — as others have before us 
using vastly different approaches — that the proton fraction in the inner crust is very small 
indeed. However, at this moment it is unclear whether exotic pasta structure can develop in 
such proton-poor environments. In order to answer this question one would have to identify 
physical observables that are sensitive to the formation of the pasta. This task, however, 
has so far proven elusive. In an effort to answer these open questions we plan to perform 
Molecular Dynamics simulations with a very large number of particles. By doing so we will 
be able to calculate both static and dynamic observables. We trust that hidden among these 
host of observables will be at least one that will respond dramatically to the formation of 
the nuclear pasta. 
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Appendix: The Ewald Construction 

In this appendix we describe the derivation of the Ewald formula used to compute the 
Coulomb potential of a system of Z protons immersed in a neutralizing and uniform electron 
background. Although such a derivation is now part of the standard repertoire of Computer 
Simulation textbooks (see for example Refs. pH - [26] ) we present here a detailed account in 
an attempt to keep the manuscript self-contained. 

We consider a system of Z point protons confined to a simulation cube of (finite) volume 
V = L^ that is immersed in a neutralizing electron background of uniform density Pe = 
—eZ/V. The Coulomb interaction for such a system is given by 

Vcoui{ri,...,rz) = - p{v)^{v)dh = -eY,^(v,) - y $(r)dV , (Al) 

Jv j^]^ Jv 

where $(r) is the electrostatic potential that one must compute as a solution of Poisson's 
equation. In an effort to minimize finite-size effects we employ periodic boundary con- 
ditions. That is, the (finite) simulation box is exactly replicated an infinite number of 
times — resulting in a perfect periodic tiling of all space. It is this resulting periodic charge 
density that must be used as the source of the electrostatic potential. In this periodic ap- 
proach, each proton within the simulation volume V interacts not only with the other Z—1 
protons contained in the finite simulation box but, in addition, with all the (infinitely many) 
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"images" residing in the periodic boxes. Clearly, such a prescription appears highly imprac- 
tical. In the case of short-range interactions, such as the strong nucleon-nucleon interaction, 
periodic boundary conditions pose no additional computational burden as each individual 
proton in the simulation box interacts with only one — the closest image — of the remaining 
Z—1 protons. This is because the next closest image must be at least half a box length {L/2) 
away, thereby suppressing the contribution from all remaining images provided, of course, 
that the range of the interaction is significantly smaller than L/2. However, the long-range 
Coulomb interaction has no intrinsic scale so one must cope with the challenge of dealing 
with all periodic images. Therein lies the power of the Ewald suniniation. 



To compute the electrostatic potential that enters in Eq. (A.l) one must solve Poisson's 
equation using the periodic charge density (including image charges) as the source term. 
That is, 

V2$(r) = -47rp(r) , (A.2) 

where 

Z 



P(r) = e E E n^ - (^^ - ^^) - ^^ • (^-3) 



=1 



V 



Here n = {nx,ny,nz) is a translational vector consisting of three arbitrary integers. To 
improve the convergence properties of the (infinite) sum required to compute the Coulomb 
interaction, Ewald introduced a distribution of smeared charges of magnitude — e centered 
around the (point) protons. Indeed, Ewald's main "trick" consists in adding and subtracting 
the following charge density: 

z 
Pair) = -e 5^ 5^ 5, (r - (r, - nL)) , (A.4) 

i=l n 

where Sa{r) represents a diffuse charge distribution (of unit charge) and smearing parameter 
a. As it is often done, we assume here a gaussian charge distribution of the form 

'^-(^) = , ,^3/2 exp(-rVa2) . (A.5) 

By adding and subtracting the above smeared charge distribution one may decompose the 
total charge density in two components; a short-range one and a long-range one. That is, 

p(r) = &r(r) + ^ir(r) , where (A. 6a) 

z 
^sr(r) = e5^5^[5(r-(r,-nL))-(5,(r-(r,-nL))] , (A.6b) 

i=l n 
Z ^ 

qM = (^^^^a[r-{ri- nL) j - e- . (A.6c) 



1. Short-range electrostatic potential 

The short-range contribution to the electrostatic potential may be directly computed 
using Poisson's equation. In particular, the electrostatic potential due to the gaussian charge 
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distribution given in Eq. (A. 5) is given by 



$a(r) 



erf(r/a) 



(A.7) 



where erf (x) is the error function. As expected, the electrostatic potential generated by such 
a smeared charge distribution modifies the Coulomb potential at short distances (r < a) but 
then "heels" back rapidly to the standard 1/r form as r ^ a. This results in an effective 
screening of the bare proton charge for distances larger than the smearing parameter. Ef- 
fectively, the long-range electrostatic potential gets replaced by a short-range one. That 
is, 



<l>(r) 



1 



-^ $sr(r) 



1 erf(r/a) _ erfc(r/a) 



(A.8) 



ly iyi itn /y 

where erfc(x) = 1 — erf(a;) is the complement of the error function. Using the principle of 
superposition, the short-range component of the electrostatic potential is given by 



$sr(r; a) 



EE 

i=l n 



erfc 



[y — Yi + nL)/a 



\y -Yi + nL\ 



(A.9) 



2. Long-range electrostatic potential 



The remaining — long-range contribution — will be computed using Fourier-transform tech- 



niques. Given that both the long-range component of the charge density [^ir(i') in Eq. (A. 6c) 
and the resulting electrostatic potential are periodic functions of r, they can be expressed 
in terms of a Fourier series as follows: 



$ir(r) = yYl ^^'^^'> exp(^k ■ r 

k 

ar(r) = yYl ^^'^^'> 6^P(^k ■ r) 



(A.lOa) 
(A.lOb) 



where the allowed momenta k are quantized in units of 2tt/L. In terms of the Fourier 
representations, Poisson's equation may be solved by inspection. That is. 



$lr(k) 



An 

^^lr(k) 



(A.ll) 



where the charge form factor f?ir(k) is given by 



.„(k)=/<;Vexp(-*.r) 

Jv 



' XI X^ ^"^^ - rj + nL) - e 



i=l n 



V 



(A.12) 



Note that as it stands, $ir(k) diverges at k = due to the long-range nature of the Coulomb 
interaction. As we will show below, the presence of a neutralizing electron background 
leads to ^ir(k = 0) = (as a long wavelength probe can only resolve the total charge of the 
system) — effectively removing this zero-momentum singularity. 

One of the great simplifications that emerges from the periodicity of the charge distribu- 



tion is that the integral over the finite simulation volume V in Eq. ( A.12 ) can be transformed 
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into an integral over all of space (Q). This fact may now be used to obtain the following 
simple and illuminating form for £iir(k): 

z „ 

^i,(k) = e^ / rfVexpHk-r)5,(r-r,)-eZ,5k,o = i^ch(k)5,(k)(l-Vo) , (A.13) 

where -Fch(k) is the charge form factor of the (point) proton distribution and ^^(k) is the 
Fourier transform of the smearing function. That is, 



Fch(k) = e ^ exp(-ik ■ r^) 



(A.14a) 



6a{k) = f d^rexp{-ik ■ r)(5„(r) = expi-k'^a^/i) . (A.14b) 

Jq 

Three features were instrumental in obtaining a form for ^ir(k) that makes the formalism 
amenable to numerical computations: (a) the periodicity of the charge distribution, (b) the 
diffuseness of the single-proton density, and (c) the charge neutrality of the system. First, 
the periodicity of the system enabled one to transform an integral over the finite simulation 
volume V over an integral over all of space, thereby simplifying the calculation of (5a (k). 
Second, the diffuseness of the single-proton density suppresses momentum components that 
are much larger than a~^, leading to a rapid convergence of the momentum sum required to 



compute $ir(r) [see Eq. (A. 15)]. Finally, the neutralizing electron background removes the 
zero-momentum component of f?ir(k), rendering the electrostatic potential finite at k = 0. 

Using Eqs. (A. 10a), (A. 11), (A.13), and (A. 14b), we obtain the following (numerically 
suitable) form for the long-range component of the electrostatic potential: 



$ir(r; a) 



Air 

v 



E^^h(k)- 



kT^O 



-fc2a2/4 



Ali.r 



(A.15) 



Having obtained the electrostatic potential using Ewald's construction, we are now in 
a position to compute the resulting Coulomb energy. As was done with the electrostatic 
potential, we separate the Coulomb energy into a short-range and a long-range component. 



3. Coulomb-Ewald potential 

The short-range component of the Coulomb potential is obtained by inserting the short- 



range electrostatic potential $81(1") given in Eq. (A.9) into Eq. (A.l). We obtain 



y3r(ri, ...,rz;a) = -eJ2 '^sr(r.) - ^ / $sr(r)t/V 

1 V z . r 

= - ^ Wsr(ri - Fj + nL) - ^ 5Z / Vsr{r - Fi + nL) d^r , (A. 16) 

where we have introduced the modified — short-range — two-body "Coulomb" interaction as 

follows: 

2 erfc(r/a) 



v^Ar = e 



(A.17) 
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The electronic contribution [second term in Eq. (A. 16)] may be evaluated by using (as we 
did earlier) the periodicity of the problem to re-write the integral over the finite simulation 
volume as an integral over all of space. That is, 



V / i;sr(r - Fi + nL) dV = V / 
,„ Jv , Jn 



■Usrfr - Yi)d r . 



(A.18) 



Given that the integral is now over all space, one then can shift the variable of integration to 
remove all dependence on the proton coordinates r^ to obtain a (constant) term proportional 
to Z: 



2=1 



v.rir - r,) d r = Z 



/ ^^sr( 

Jn 



r) d^r = Ze^Tra^ . 



(A.19) 



As it stands, the sum over proton coordinates in Eq. (A. 16)] diverges because of spurious 



self-interactions (i.e., the Z terms with i = j and n = 0). To render this sum finite we must 
remove all self-interactions among the point protons. That is. 



K 



self 



-Z lim — 

2 r^o r 



Ze^ 
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2 r-5>0 
Z 



erfc(r/a) erf(r/a) 
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i=l 



Ze' 



a\/Tr 



(A.20) 



By subtracting this self-interaction term from Eq. (A. 16), we obtain the short-range contri- 
bution to the electrostatic Coulomb energy: 



, Iv^, , ,, Ze^ Tra^Ze)^ 
Vsr{ru...,rz;a) = ^ 2^ v^^irij + nL) - —-= - — -— , 



«J,n 



^^^Zef 
2 V 



(A.21) 



where rjj = rj — r^ and the prime is used to indicate that self-interaction terms (namely, the 
Z terms with i=j and n = 0) should all be omitted from the sum. 

As in the case of the short-range interaction, the long-range component of the Coulomb 



potential is obtained by inserting the long-range electrostatic potential $11(1") [Eq. (A. 15) 
into Eq. (lA.ll). We obtain 



V^r(ri, . . . , r^; a) = -e ^ <l>i,(ri) - ^ / $ir(r)rfV 



2V .,y 



Y, |F,h(k)p ^^^ - — $i.(k = 0) . 



V 



k^O 



A;2 



W 



(A.22) 



The first term in the above equation is obtained by directly substituting $ir(r) into the 
proton sum and then using the definition of the (point) charge form factor introduced in 
Eq. (A. 14a). The second term is proportional to the zero- momentum component of the 



potential and vanishes because of the overall charge neutrality of the system [see Eq. (A. 13 )] 
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Collecting both — short- and long-range — contributions we arrive at a form of the Coulomb 
potential that is both well defined and rapidly convergent: 

Vboui(ri, ■■■,rz) = yYl l^=h(k)P ^, + - J2'vsM, + nL) 

a,/^ 2 V ■ ^^■^'^' 



Before leaving the derivation of the Ewald method and proceed to test its accuracy and 
convergence properties, we present the explicit form of the Coulomb potential that will be 
used in our simulations. To do so we take advantage of the fact that the long-range Coulomb 
potential is an interaction with no intrinsic scale. To start, we use the box-size L to define 
dimensionless coordinates (s) and momenta (1) as follows: 

r = Ls , (A.24a) 

27r 
k = — 1 . (A.24b) 



By introducing these simple definitions into Eq. (A. 23), we isolate the dimensionfuU pa- 
rameter fo = e^/L from the dimensionless many-body Coulomb potential. That is, 

Vboui(ri, ...,rz)= i — ) f/coui(si, . • . , Sz) , (A. 25) 



where 



...,rz) = 


- -y f^Coui(si,. 


^E^=. 

¥0 


.(i)'MO + 





(A.26) 



Here the long-range interaction is written in terms of the (point-proton) charge form factor 
Fch(l) and the modified Coulomb propagator w{l) 

z 
Fch(l) = ^ exp(-27ril ■ s^) , (A.27a) 



exp( — TT^Sg/^) 
7/2 



w{l) = ^^ _ ° > (with So = a/L) , (A.27b) 



the short-range interaction is given by 



.„(s) = ?':f!!<^/5?) . (A,28) 



and U'q is a constant: 



C/; = -^-^. (A.29) 

Evidently, the Coulomb potential must be independent of the auxiliary smearing param- 
eter sq. This is convenient as one may use this independence to test the reliability of the 
method. Once such a test has been performed, one can use this freedom by selecting a value 
of So that will expedite the computation without sacrificing accuracy. In particular, in this 
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contribution we will employ a smearing parameter that is significantly smaller than the box 
length (i.e., Sq ?« 0.1). This ensures that the mininium image convention may be used to 
compute the short-range part of the Coulomb potential. To compute the long-range part of 
the potential one must evaluate a sum over momentum modes. For a small value of Sq, as 
the one adopted here, a very large number of modes must be included in order to achieve 
convergence. As a result, this part of the computation can be quite expensive. Thus, one 
often resorts to methods that handle the Fourier sum more efficiently, such as particle-niesh 
approaches [25j. In this contribution we implement a simple "poor's-man" approach that 
requires pre-computing the potential on a relatively fine mesh. During the simulation, com- 
puting the long-range part of the potential is then reduced to a simple interpolation scheme. 



To do so we write the square of the of the charge form-factor appearing in Eq. (A. 26) as 

follows: 

z z z 

|Feh(l)P = ^^exp(-27ril ■ Sij) = Z + ^exp(-27rzl ■ s^j) . (A.30) 

In this way the long-range component of the Coulomb potentials reduces — as in the case of 
the short-range component — to a sum of two-body terms (plus a constant). That is, 

f/ir(si, ...,sz) = ^Yl l^ch(l)|' w(0 = ^ $^Mir(si,) + |nir(0) , (A.31) 

where the (long-range) two-body interaction has been defined by 

Mj^(s) = Y^ w{l) exp(-27rzl • s) . (A.32) 

¥0 

Finally, the dimensionless Coulomb interaction to be used in the Monte Carlo simulations 
is given by the following expression: 

f^Coul(Sl, . . . , S^) = - ^ I^Mlr(Sij) + Msr(Sij) + Uq ■ (A. 33) 



¥i 



C/„ = |.„(0) - ^ - !^ . (A.34) 

2 V^-^o 2 



where 



The Ewald construction is a powerful method that may be tested in a variety of ways. 
First, one may compare against well known results for the Coulomb energy of simple crys- 
tals. Second, whereas the Ewald construction introduces a dependence on the smearing 
parameter a, the "bare" Coulomb potential should be insensitive to it — as long as a is nei- 
ther too big nor too small. We have implemented such ideas in Fig. [7] where the Coulomb 
energy of a system consisting of Z protons localized to the sites of a simple cubic lattice 
and immersed in a uniform electron background has been computed. The figure displays 
the competition between the long-range (repulsive) and the short-range (attractive) contri- 
butions to the Coulomb energy as a function of the smearing parameter a and number of 
protons. The energy is given in units of e^/aiatt, where aiatt = Pz i^ ^^^ lattice constant and 
Pz is the proton density. We obtain for all values of a and all system sizes an energy per 
particle of —2.837297479/2, where 2.837297479 is the Madelung constant for a simple-cubic 
crystal [21] • Remarkably, the Ewald method can already reproduce the exact answer (up to 
nine significant figures!) for a system with two protons per side. 
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FIG. 7: (Color online) Coulomb energy per particle together with its long- and short-range con- 
tributions for a simple cubic lattice of Z protons immersed in a neutralizing and uniform electron 
background. The Coulomb energy is given in units of e^/aiatt (oiatt is the lattice constant) and the 
dimensionless Madelung constant for a simple cubic lattice is equal to 2.837297479 [31]. 
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